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QUADRATIC OPTIMIZATION IN THE PROBLEMS OF ACTIVE CONTROL OF 

SOUND* 


J. LONCARICt AND S. V. TSYNKOV* 

Abstract. We analyze the problem of suppressing the unwanted component of a time-harmonic acoustic 
field (noise) on a predetermined region of interest. The suppression is rendered by active means, i.e., 
by introducing the additional acoustic sources called controls that generate the appropriate anti-sound. 
Previously, we have obtained general solutions for active controls in both continuous and discrete formulations 
of the problem. We have also obtained optimal solutions that minimize the overall absolute acoustic source 
strength of active control sources. These optimal solutions happen to be particular layers of monopoles on 
the perimeter of the protected region. Mathematically, minimization of acoustic source strength is equivalent 
to minimization in the sense of L\. 

By contrast, in the current paper we formulate and study optimization problems that involve quadratic 
functions of merit. Specifically, we minimize the L 2 norm of the control sources, and we consider both 
the unconstrained and constrained minimization. The unconstrained L 2 minimization is certainly the eas- 
iest problem to address numerically. On the other hand, the constrained approach allows one to analyze 
sophisticated geometries. In a special case, we cam compare our finite-difference optimal solutions to the 
continuous optimal solutions obtained previously using a semi-analytic technique. We also show that the 
optima obtained in the sense of L2 differ drastically from those obtained in the sense of L\. 

Key words, noise cancellation, active control sources, volumetric and surface controls, general solution, 
monopoles and dipoles, radiation of waves, complex-valued quantities, L2-minimization, overdetermined 
systems, least squares, unconstrained minimization, constrained minimization 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. In the simplest possible formulation, the problem of active control of sound is posed 
as follows. Let H C be a given domain (bounded or unbounded), and T be its boundary: T = <9fl, where 
the dimension of the space n is either 2 or 3. Both on Q and on its complement Q\ = R n \H we consider the 
time-harmonic acoustic field u = u(®), iGl", governed by the non- homogeneous Helmholtz equation: 

Lu = A u + k 2 u = /. (1.1) 
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Equation (1.1) is subject to the Sommerfeld radiation boundary conditions at infinity, which for n ~ 2 are 
set as 


(z) = o(|z| 1/2 ) , +iku(x) = o(|*| 1/2 ), as |*| 


(1.2a) 


and for n = 3 as 


u{x) = O (|*| *), -finr +iku(x) =o(|*| *) as |z| — ► oo. 


(1.2b) 


The Sommerfeld boundary conditions specify the direction of wave propagation, and distinguish between 
the incoming and outgoing waves at infinity by prescribing the outgoing direction only; they guarantee the 
unique solvability of the Helmholtz equation (1.1) for any compactly supported right-hand side / = /(sc). It 
is important to mention that as we are dealing hereafter with the traveling waves (radiation of sound toward 
infinity), all the resulting solutions will necessarily be complex- valued, otherwise it is impossible to account 
for the key phenomenon of variation of phase with the change of spatial location. 

The source terms / = f(x) in equation (1.1) can be located on both fi and its complement Cl\ = M n \ft; 
to emphasize the distinction, we denote 


/ = / + + / , supp/+ C Cl, supp / Cfii- 


(1.3) 


Accordingly, the overall acoustic field u = u(x) can be represented as a sum of the two components: 


U = U + -1- u , 


(1.4) 


where is driven by the interior sources Z -1 ", and u is driven by the exterior sources / w.r.t. Cl: 


Lu+ = /+, (1.5a) 

Lu~ = /“. (1.5b) 


Note, both u+ = u + (x) and u~ — u~(x) are defined on the entire R n , the superscripts and ” refer 
to the sources that drive each of the field components rather than to the domains of these components. The 
setup described above is schematically shown in Figure 1.1 for the case of a bounded domain Cl. 
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Fig. 1.1. Geometric setup. 


Hereafter, we will call the component u + of (1.4), (1.5a) sound , or “friendly” part of the total acoustic 
field; the component u~ of (1.4), (1.5b) will accordingly be called noise, or “adverse” part of the total 
acoustic field. In the formulation that we axe presenting, fi will be a predetermined region of space to be 
protected from noise. This means that we would like to eliminate the noise component of u(x) inside fi, while 
leaving the sound component there unaltered. In the mathematical fr amework that we have adopted, the 
component u~ of the total acoustic field, i.e., the response to the adverse sources f~ [see (1.3), (1.4), (1.5)], 
will have to be canceled out on fi, whereas the component u + , i.e., the response to the friendly sources /+, 
will have to be left unaffected on fi. A physically more involved but conceptually easy to understand example 
that can be given to illustrate the foregoing idea, is that inside the passenger compartment of an aircraft 
we would like to eliminate the noise coming from the propulsion system located outside the fuselage, while 
not interfering with the ability of the passengers to listen to the inflight entertainment programs or simply 
converse. Another good example is found in medicine, where high levels of periodic noise are produced by 
resonance coils in magnetic resonance imaging (MRI) machines. 

The concept of active noise control that we will be discussing implies that the component u~ is to be 
suppressed on fi by introducing additional sources of sound g = g(x) exterior with respect to fi, supp g C fii, 
so that the total acoustic field u = u(x) be now governed by the equation [cf. formulae (1.1), (1.3)]: 

Lu- / + + /“ + g, (1.6) 


and coincide with only the friendly component on the domain fi: 


= 

xen xeQ’ 


(1.7) 


The new sources g = g(x) of (1.6), see Figure 1.1, will hereafter be referred to as the control sources or 
simply controls. An obvious solution for these control sources is g = —f~. This solution, however, is clearly 
sub-optimal because on one hand, it requires an explicit and detailed knowledge of the structure and location 
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of the sources /“, which is, in fact, superfluous, see [4]. On the other hand, its implementation in many 
cases, like in the previously mentioned example with an airplane, may not be feasible. Fortunately, there 
are other solutions of the foregoing noise control problem (see Section 2, as well as our previous work [4] for 
detail), and some of them may be preferable from both the theoretical and practical standpoint. 

To conclude the introduction, let us only mention that the area of active control of sound has a rich 
history of development, both as a chapter of theoretical acoustics, and in the perspective of many different 
applications. It is impossible to adequately overview this extensive area in the framework of a focused 
research publication. As such, we simply refer the reader to the monographs [1,2,7] that, among other 
things, contain a detailed survey of the literature. Potential applications for the active techniques of noise 
control range from the aircraft industry to manufacturing industry to ground and air transportation to the 
military to consumer products and other fields, including even such highly specialized and narrow areas as 
acoustic measurements in the wind tunnels. It is generally known that active techniques are more efficient 
for lower frequencies, and they are usually expected to complement passive strategies (sound insulation, 
barriers, etc.) that are more efficient for higher frequencies, because the rate of sound dissipation due to 
viscosity of the medium and heat transfer is proportional to the square of the frequency [3]. 

Let us also note that in the current paper we focus on the case of the standard constant-coefficient 
Helmholtz equation (1.1), which governs the acoustic field throughout the entire space R n . This allows 
us to make the forthcoming analysis most straightforward. However, one can as well consider other, more 
complex, cases that involve variable coefficients, different types of far-field behavior, discontinuities in the 
material properties, and maybe even nonlinearities in the governing equations over some regions. Approaches 
to obtaining solutions for active controls in these cases are based on the theory of generalized Calderon’s 
potentials and boundary projections, and can be found in our previous paper [4] and in the monograph by 
Ryaben’kii [12, Part VIII]. 

The material in the rest of the paper is organized as follows. In Section 2, we introduce and discuss general 
solutions for controls in the continuous and discrete framework. Section 3 is devoted to the formulation and 
solution of the quadratic optimization problems for the control sources (unconstrained and constrained L 2 
optimization). For reference purposes we also briefly mention our previous results on the optimization in 
the sense of Lx . Finally, Section 4 provides a summary and outlines a perspective for the future work. 

2. General Solutions for Control Sources. 

2*1. Continuous Formulation of the Problem. A general solution for the volumetric continuous 
control sources g = g(x) is given by the following formula (Hi = M n \fl): 

g(x) = -Lw\ x€Ui , (2.1) 


where w = tc(x), x £ fix, is a special auxiliary function-parameter that parameterizes the family of controls 
(2.1). The function w(x) must satisfy the Sommerfeld boundary conditions (1.2a) or (1.2b) at infinity, and at 
the interface T, the function w and its normal derivative have to coincide with the corresponding quantities 
that pertain to the total acoustic field u given by formula (1.4): 



dw du 
dn r dn T 


( 2 . 2 ) 
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Other than that, the function w(x) used in (2.1) is arbitrary, and consequently formula (2.1) defines a large 
family of control sources, which provides ample room for optimization. The justification for formula (2.1) 
as the general solution for controls can be found in [4]. In our recent paper [5] we also emphasize that the 
controls 

g(x) = J g(y)S(x - y)dy = g*S 

given by (2.1) are actually volumetric control sources of the monopole type with regular density g 6 L\ ^(E n ) 
[assuming that w(x) was chosen sufficiently smooth to guarantee local absolute integrability of g(x)]. 

The control sources (2.1) possess several important properties. First of all, we see that to obtain these 
controls one needs no knowledge of the actual exterior sources of noise In other words, neither their 
location, nor structure, nor strength are required. All one needs to know is u and ^ on the perimeter T 
of the protected region ft. In a practical setting, u| r and |^| r can be interpreted as measurable quantities 
that are supplied to the control system as the input data. Let us emphasize that the quantities to be 
measured refer to the overall acoustic field u rather than only to its unwanted component u~ , see formula 

(2.2) . At the same time, the analysis of [4] shows that the application of the controls (2.1) will result in the 
cancellation of only the adverse noise u~ on the protected domain ft, whereas the friendly sound field u+ 
will be left unaffected. In other words, the controls (2.1) are insensitive to the interior sound u + , whatever 
it might be, and are built so that to suppress only the exterior noise u~ on ft. This capability is extremely 
important because in many applications the overall acoustic field always contains a component that needs to 
be suppressed along with the part that needs to be left intact. Let us also note that a more general analysis 
of [4] based on Calderon’s potentials and boundary projections yields the same formula for controls (2.1), 

(2.2) for the cases that may involve variations in material properties and alternative types of the far-field 
behavior. Of course, the operator L will be a new variable-coefficient operator, and the function-parameter 
u>(a:) will have to satisfy new far-field boundary conditions instead of the Sommerfeld boundary conditions. 

Along with the volumetric controls (2.1), one can also consider surface controls, i.e., the control sources 
that are concentrated only on the interface F. A general solution for the surface controls is given by 

— ^ (2 - 3) 

where w = w(x), as before, denotes the auxiliary function-parameter. In contradistinction to the previous 
case, now it has to satisfy the homogeneous Helmholtz equation on the complementary domain: Lw = 0 for 
x € fti, and the Sommerfeld boundary condition (1.2a) or (1.2b) at infinity, but at the interface T it may be 
arbitrary, i.e., it does not have to meet boundary conditions (2.2). The corresponding discontinuities that 
are denoted by expressions in rectangular brackets in formula (2.3) drive the surface control sources. The 
first term on the right-hand side of (2.3) represents the density of a single-layer potential, which is a layer 
of monopoles on the interface T, and the second term on the right-hand side of (2.3) represents the density 
of a double-layer potential, which is a layer of dipoles on the interface F. A detailed justification of formula 

(2.3) as general solution for surface controls can be found in [15], see also [5]. The fundamental properties of 
the surface controls (2.3) are the same as those of the volumetric controls (2.1) — they are also insensitive 
to the interior sound u+(x), and do not require any knowledge of the actual sources of noise /". 

In the family of surface controls (2.3) we identify two important particular cases. First, the cancellation 
of u~(x), x € ft, can be achieved by using surface monopoles only, i.e., by employing only a single-layer 
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potential as the annihilating signal (anti-sound). To do that, we need to find w(x ), x E fix, such that there 
will be no discontinuity on T between u(ar) and w(x), i.e., in the function itself, and the discontinuity may 
only “reside” in the normal derivative [see formula (2.3)]. This w(x) will obviously be a solution of the 
following external Dirichlet problem: 


Lw = 0, x E Hi, 


w 


lr 



(2.4) 


subject to the appropriate Sommerfeld boundary condition (1.2a) or (1.2b). Problem (2.4) is always uniquely 
solvable on = R n \H. Second, one can employ only the double-layer potential to cancel out u~(x), x E fi, 
i.e., use only surface dipoles as the control sources. In this case, the function w(x), x E fli, has to be chosen 
such that the discontinuity on T be only in the function itself, i.e., between the actual values of u{x) and 
w(x ), and not between the normal derivatives. This w(x) should then solve the following external Neumann 
problem: 


Lw = 0, x E fix, 


dw 

Ihi 


du 

dn 


(2.5) 


again, subject to the appropriate Sommerfeld condition at infinity, (1.2a) or (1.2b); the latter guarantees 
the solvability of (2.5). We therefore see that surface control sources (2.3) are given by combinations of the 
monopole and dipole layers, with the two “extreme” cases corresponding to either only monopoles, see (2.4), 
or only dipoles, see (2.5). 

Altogether, we have now introduced active controls of two different types on the surface, but only one 
type of the volumetric controls — monopoles, see formulae (2.1), (2.2). This is not accidental. Let us note 
that from the standpoint of physics and engineering, the monopole and dipole sources provide different types 
of excitation to the surrounding sound-conducting medium. A point monopole source can be interpreted as 
a vanishingly small pulsating sphere that radiates acoustic waves symmetrically in all directions, whereas a 
dipole source resembles a small oscillating membrane that has a particular directivity of radiation. Moreover, 
in the genuine time-dependent acoustic context one can show that monopole sources are those that alter the 
balance of mass in the system, they are scalar in nature and reside on the right-hand side of the continuity 
equation, whereas dipole sources alter the balance of force, they are vectors and reside on the right-hand 
side of the momentum equation, see our recent work [5] for detail. This distinction basically warrants a 
separate consideration of the monopole and dipole type sources as far as the point-wise or surface excitation 
may be concerned. As, however, has been shown in [5], in the framework of time-harmonic volumetric 
excitation (the case studied hereafter) a separate consideration of dipole fields appears superfluous. In fact, 
any volumetric distribution of dipoles can, under the assumption of sufficient regularity, be recast in the form 
of an equivalent volumetric distribution of monopoles. In so doing, the dipole sources enter the right-hand 
side of the Helmholtz equation (1.1) through a divergence operator, whereas monopoles enter this right- 
hand side directly (up to a multiplicative constant). 1 We refer the reader to our paper [5], as well as to 
the monograph [7], for further detail. In Section 3, we will study the volumetric monopole controls in the 


1 In the corresponding analysis, we interpret the field variable it(x) as acoustic pressure, which is a common strategy in the 
literature. 
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context of L 2 optimization; for comparison, we also provide there the results of the Li optimization from [5] 
that involve both the volumetric and surface monopole control sources. 

Let us also note that in practice it may often be convenient to use the so-called artificial boundary 
conditions (ABCs), see [14], as a part of the definition of the auxiliary function w(x). Assume that there is 
a larger domain that fully contains fl and require, in addition, that Lw = 0 outside this larger domain. This 
requirement is always met in the case of the surface controls (2.3); and in the case of the volumetric controls 
it implies that the resulting control sources will be compactly supported between T and the outer boundary 
of the aforementioned larger region, see formula (2.1). For many applications this is desirable. Moreover, 
from the standpoint of computing this is clearly the only feasible way to obtain a finite discretization (see 
Section 2.2). It is known that the homogeneous equation Lw — 0 outside a given region, along with the 
Sommerfeld boundary conditions at infinity, can be equivalently replaced by special ABCs at the boundary 
of this region. General approaches to building the ABCs for a variety ol different formulations are discussed 
in the review paper [14]. For the specific case of a homogeneous Helmholtz equation outside a sphere of 
radius R in 3D, the ABCs were obtained in [5] using the separation of variables in spherical coordinates and 
mode selection that would guarantee that the boundary conditions at infinity are satisfied: 


dwim 

dp 


^{p- 1/2 Hjl\ /2 (kp)] „ 


/-f 1/2 ‘ 


I p=R 


p- 1/2 H,% 2 (kp) 


-Wlrn(p) 


I p=R 


( 2 . 6 ) 


In formula (2.6), p is the spherical radius, w lm are the Fourier coefficients of w(x) with respect to spherical 
functions Y t m , l = 0,1,2,..., m = 0, ±1, . . . , ±1, and h\+ 1/2 are Hankel’s functions of the second kind; 
equalities (2.6) have to be enforced for all the appropriate l and m. Similarly, for the homogeneous Helmholtz 
equation outside a disk of radius R in 2D, the ABCs obtained in [5] read: 


dwi 

dp 


4,H\ 2 \kp) „ 




h[ 2) 


(kp) 


wi(p) 


(2.7) 


p=R 


where p is the polar radius, and wi are the Fourier coefficients of w(x) with respect to the complex exponents 
e~ il0 , l = 0, ±1, ±2, . . again, equalities (2.7) have to be enforced for all L 


2.2. Discrete Formulation of the Problem. The continuous analysis tools employed for obtaining 
the control sources of the previous Section 2.1 are obviously deficient from the standpoint of applications. 
Indeed, any practical design of a noise control system can only be composed of a finite number of elements 
(sensors for measuring the field and actuators for creating the appropriate excitation, i.e., anti-sound). 
Therefore, it is natural to discretize the problem on the grid and obtain the control sources in the discrete 
framework so that the locations of the sensors and actuators can be associated with the grid nodes. For 
details regarding the discrete formulation of the problem we refer the reader to the monograph [12, Part 
VIII], as well as to the papers [17, 18]; a brief account can also be found in [5, 15], and below we summarize 
the results. Note that our discrete analysis is not limited to any specific type of the grid. In particular, 
no adaptation or grid fitting to either the shape of the protected region ft or that of the external artificial 
boundary, is generally required. However, for the purpose of illustrating the concepts discussed hereafter, we 
will use a two-dimensional example that involves a polar grid. The use of the polar grid greatly facilitates 
setting the discrete ABCs at the circular outer boundary of radius R . Moreover, the same two-dimensional 
polar example is analyzed later in Section 3 in the context of Li optimization. 
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Let us denote the aforementioned polar grid N; it spans both ft and ft\. Of course, the grid does not 
extend all the way to infinity, it is rather truncated by the external artificial boundary in the shape of a large 
circle of radius R . This, in particular, implies that the discrete control sources that we obtain will always 
be compactly supported (see the discussion in the end of Section 2.1). Assume that the grid has J cells in 
the radial direction with the nodes pj = jAp, j = 0, . . . , J, so that po = 0 and pj = R\ and L cells in the 
circumferential direction with the nodes 6 S = sA0, s = 0, . . . , X, so that 6q = 0 and Ql = 27i\ For simplicity, 
it is convenient to think that the grid sizes A p — R/J and AO = 2n jL are constant; in applications, however, 
the grid in the radial direction may be stretched. 

Let now be a representation of the acoustic field on the grid, and L ^ be a finite-difference ap- 
proximation of the differential operator X of (1.1). To accurately define the approximation, we will need to 
introduce another grid M along with the previously defined N. On the grid M, we will consider the residuals 
of the operator lS h \ and subsequently the right-hand sides to the corresponding discrete inhomogeneous 
equation. We will use the notations n and m for the individual nodes of the grids N and M, respectively, 
and the notation Km for the stencil of the discrete operator X^ centered at a given node m € M, so that 


L (h) u {h] 


m 


E a u w 
u mn u ’n > 

n€N m 


( 2 . 8 ) 


where a nm are the coefficients associated with particular nodes of the stencil. Generally, there are no 
limitations to the type of the discrete operator that one may use. We only require that the difference 
operator X^ of (2.8) approximate the differential operator X of (1.1) with the accuracy sufficient for a 
particular application. For the specific example that we are analyzing, we will consider a conventional 
second-order central- difference approximation, so that the grids N and M actually coincide: n — (s, j) and 
m = ( s,j ) (except that M is smaller, it does not contain the outermost row p — pj = R), and formula (2.8) 
becomes: 


2,W u (fc)l . = IJ_ 

ls >J pj Ap 


Pj+i 


,(h) 

Vi+i 


i (h > 


Ap 


-pj 


u {h) - u (h) 
U sJ U 8,j - 1 




A p 


) 


+ 


i {h) 


ZM 




A 9 2 


i {h) 

+ k 2 u[ h) j. 


(2.9) 


Next, we introduce the following subsets of the grids M and N, which will allow us to accurately dis- 
tinguish between the interior and exterior domains, interior and exterior sources, and interior and exterior 
solutions on the discrete level: 


M + =Mnn, MT =M\M + =Mnfii, 

N+ = U N m , N- = U N m , 

mEM 4 mEM - 

7 = nr, 7 + = r nfi, y- = N^nftx. 


( 2 . 10 ) 


We emphfisize, that the grid M that pertains to the residuals of the finite-difference operator L ^ is parti- 
tioned into M + and M - directly, i.e., following the geometry of ft and ft\. In contradistinction to that, the 
grid N is not partitioned directly, we rather consider the collection of all nodes of N swept by the stencil 
N m when its center 2 belongs to fl, and call this sub-grid N + , see (2.10). Obviously, some of the nodes of 

2 These definitions require that the grid be sufficiently large to ensure that the sets M+ and M - are not empty; each must 
include at least one layer of nodes along the interface, as in formula (3.11), see Section 3.2. 
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Fig. 2 . 1 . Schematic geometry of the domains, the stencil , and the grid boundary 7 = 7+U7 in polar coordinates: Hollow 
bullets denote 7+, filled bullets denote 7". 


N + obtained by this approach happen to be outside ft, i.e., in fti, and these nodes are called 7“. The sets 
N" and 7+ are defined similarly starting from M - . The key idea is that whereas the grids M + and M~ do 
not overlap, the grids N+ and N - do overlap, and their overlap is denoted 7; obviously, 7 = 7+ U 7". The 
subset of grid nodes 7 is called the grid boundary, it is a fringe of nodes that is located near the continuous 
boundary T and in some sense straddles it. The specific structure of 7 clearly depends on the construction of 
the operator Z^ of (2.8) and the stencil N m . For example, for the polar second-order Laplacian (2.9), the 
grid boundary 7 will be a two-layer fringe of grid nodes located near T, as shown schematically in Figure 2.1. 
Further specifics on the construction of grid boundaries can be found in the monograph [12]. 

The discrete noise control problem is formulated similarly to the continuous one, see Section 1. Let /^ + , 
m G M + , and /m m £ M~ , be the interior and exterior discrete acoustic sources, respectively. Let 
n G N, and n € N, be the corresponding solutions, i.e., Z^W^ 4- = /^ + and 

Using the same terminology as before, we will call the discrete sound and u^~ the discrete noise. 

The overall discrete acoustic field is the sum of its sound and noise components, + u^~ on 

N, and obviously satisfies the equation Z( fc )uW = /W == /^ 4 + The goal is to obtain the discrete 

control sources pW = so that the solution of the equation L^ h -u^ = + f( h )~ + <7^ be equal 

to only the sound component tt^ + on the sub-grid N + . 

A general solution for the discrete control sources g ^ = gl£^ that eliminate the unwanted noise 
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on N + is given by the following formula [cf. formula (2.1)]: 


gW = - L W W W 


m€M“ 


(2.11) 


where = Wn , n E N _ , is a special auxiliary grid function-parameter that parameterizes the family of 
controls (2.11). The requirements that this function must satisfy are, again, rather “loose,” and can 
be considered natural discrete counterparts of the corresponding requirements of the continuous function- 
parameter w(x). Namely, at the grid boundary 7 the function has to coincide with the overall acoustic 
field u M to be controlled: 


w 


(h) 

n 


= **< k) • 
n € 7 n 7 


( 2 . 12 ) 


Notice that since, e.g., for the second-order discretizations the grid boundary 7 contains two layers of nodes, 
7 + and 7", see Figure 2.1, then specifying the corresponding nodal values on 7 is in some sense equivalent 
to specifying the function and its normal derivative on T in the continuous case, see (2. 2). 3 We also note 
that for practical designs, the boundary data t4^| n€7 shall be interpreted as measurable quantities that 
provide input for the control system. In other words, we can think of a microphone at every node of 7; these 
microphones measure the characteristics of the actual acoustic field and generate the input signal | nG7 - 
The other requirement of the function w^ h \ besides the interface boundary conditions (2.12), is that 
they must satisfy the appropriate discrete ABCs at the external artificial boundary p = R, see Figure 2.1. 
The role of the discrete ABCs is the same as that of the continuous ABCs — to provide a replacement for 
the Sommerfeld radiation boundary conditions. The discrete two-dimensional ABCs are obtained in [5] by 
using the direct and inverse discrete Fourier transforms, / = —L/2 + 1, . . . , L/ 2, s = 0, . . . , L — 1: 


Wl 


L-l 


1 

E 


w 8 e 


-its AO 


8=0 


L/2 

w B = y, 

l=-L/ 2+1 


(2.13) 


and essentially approximating boundary conditions (2.7) for l = — L/2 + 1, . .. , L/2 with the second order 
of accuracy: 


u>i,J ~ m f J - 1 

Ap 


-Pi 




= 0, = 


t p h, 


( 2 ) 

Ctl 


(Ml 


H^(kp) 


OL 


2 _ 


I p=R 


Ad 2 


sin' 


IA6 

2 


(2.14) 


Note, —af of (2.14) are eigenvalues of the circumferential component of the discrete Laplacian (2.9). Relations 
(2.14) for all l = —L/2 -hi,..., L/2 can also be recast into the matrix form: 


where F and F~ l are matrices of the direct and inverse discrete Fourier transforms of (2.13), and w.j and 
w.j-i are L-dimensional vectors of components w[ h j and respectively, s = 0, 1, . . . , L - 1. 

Other than the two aforementioned requirements, i.e., the interface conditions (2.12) and the ABCs 
(2.15), the function is arbitrary and as such, parameterizes a substantial variety of discrete control 
sources, see (2.11). The latter will provide the search space for optimization in Section 3. 


r 


F = T w. 


,j- i> 


(2.15) 


w.j = F 1 diag 


Hi 


3 This statement can be given a rigorous formulation in terms of approximation, see [12] for detail. 
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It is also important to understand in what sense this discrete cancellation of noise models the continuous 
cancellation described in Section 2. This is basically the question of approximation of the continuous gener- 
alized Calderon’s potentials by their discrete counterparts. To that effect, the theory of difference potentials, 
see [12], says that under certain natural conditions, the discrete anti-sound = v n^, n £ N+ , i.e., the 
solution to = gW with given by (2.11), approximates the continuous anti-sound v = v(x ), 

x € H, i.e., the solution to Lv = g with g given by (2.1). The aforementioned natural conditions include 
first the consistency and stability of the finite-difference scheme for the Helmholtz equation. Consistency 
and stability will guarantee convergence as the grid size vanishes. In addition, the discrete boundary data 
(2*^) ^ ave to a PP rox i ma te the continuous boundary data (u, | r of (2.2) in the following 
sense. Once the continuous function u and its first-order normal derivative ^ are known at the boundary 
T, normal derivatives of higher orders can be obtained via the differential equation itself, and the near- 
boundary values u ^ | n€ can then be calculated using Taylor’s expansion; the order of accuracy of the latter 
calculation with respect to the grid size h has to be at least as high as the order of accuracy of the interior 
scheme. In this case, the quality of approximation, i.e., the rate of convergence of the discrete potential to 
the continuous one with respect to h, will be the same as prescribed by the finite-difference scheme itself. 
For the central-difference operator (2.9), this rate is 0(h 2 ). In other words, when designing an active control 
system following the finite-difference approach, one can expect to have the actual noise cancellation in the 
same approximate sense as the solution of the finite-difference equation approximates the corresponding 
solution of the original differential equation. Note, in any particular practical setting we will need to require 
sufficient wave resolution on the grid, i.e., the waves of length A = 2n/k, where k is the wavenumber in (1.1), 
will have to be well resolved by the specific discretization. 

Finally, similarly to the continuous case we can identify some particular types of the discrete control 
sources. First, let us define another subset of the grid M (more precisely, of M~): 


Mjnt = { m € M |N m n = 0} . 

Basically, M i “ t is the interior subset of M“ , such that when the center of the stencil sweeps this subset, the 
stencil itself does not touch 7 + , see Figure 2.1. In other words, M^ t is a subset of M“ such that 


u N m =N-\ 7 + . 

m GMT nt 

Next, we introduce the auxiliary function = w n 6 N~, for (2.11) as follows: 


and 


w 


W 


n£ 7 + 



n67 + ’ 


(2.16a) 


w. 


(ft) I 


= u 


(ft)l 


nG 7 — — n fnG7 — ' 

£(ftVft) =o on MT t . 


(2.16b) 


As before, we also assume that u :ih] satisfies the discrete ABCs (2.15). Definition (2.16a) means that on the 
interior part of the grid boundary 7 + we simply set equal to the given u ( h 1 : | n67+ = “i/ 1 ' | n£7+ - 

Definition (2.16b) is actually a discrete exterior boundary-value problem of the Dirichlet type. Indeed, 
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everywhere on and “outside” the exterior part of the grid boundary 7“, i.e., on N _ \7 + , the grid function 
is obtained as a solution of the homogeneous equation = 0 (enforced at the nodes Mj“ t ) 

supplemented by the boundary data on 7“: \ n€l - = u n^| n€ -> which is specified for the unknown 

function w w itself. Note, relation (2.16a) and the first relation (2.16b) together are obviously equivalent to 
(2.12). Therefore, the function w ^ defined via (2.16a), (2.16b) falls into the general class of used for 

obtaining the discrete control sources (2.11). 

Problem (2.16b) can clearly be considered a finite-difference counterpart to the continuous Dirichlet 
problem (2.4). Therefore, its is natural to call the control sources = ^monopoie obtained by formulae 
(2.11), (2.16a), (2.16b) the discrete surface monopoles . Indeed, because of the definition of w ^ given by 
(2.16a) and (2.16b), these ^monopoie ma y> generally speaking, differ from zero only on the grid set M"\M i “ t , 
which is a single “curvilinear” layer of nodes of grid M that follows the geometry of T. Accordingly, the 
output of these controls can be called the discrete single-layer potential. The discrete surface monopoles and 
discrete single-layer potential were first introduced and analyzed in our recent paper [15]. As shown in [5], 
the controls of this particular type play a key role in the context of L\ optimization, see also Section 3.1. 
Let us emphasize that unlike the continuous surface controls (2.3), which belong to a different class of 
functions rather than the volumetric sources (2.1) (singular 8 - type distributions vs. regular locally integrable 
functions), the foregoing discrete surface monopoles belong to the same original class of discrete control 
sources (2.11). Let us also note that besides the discrete surface monopoles and the corresponding single- 
layer potential, one can also define the discrete surface dipoles and, accordingly, the double-layer potential, 
see [15] for detail. 

3. Optimization of the Control Sources. Once the general solution for controls is available, in 
either continuous (2.1) or discrete (2.11) formulation, the next step is to decide what particular element of 
this large family of functions will be optimal for a specific setting. There is a multitude of possible criteria for 
optimality that one can use. In many practical problems the cancellation of noise is only approximate and as 
such, the key criterion for optimization (or sometimes, the key constraint) is the quality of this cancellation, 
i.e., the extent of noise reduction. In contradistinction to that, in this paper we are considering ideal, or 
exact, cancellation, i.e., every particular control field from either the continuous (2.1) or discrete (2.11) family 
completely eliminates the unwanted noise on the domain of interest. Consequently, the criteria for optimality 
of the controls that we can employ will not include the level of the residual noise as a part of the corresponding 
function of merit, and should rather depend only on the control sources themselves. At a later stage of the 
work we plan to look into the issues of approximate, rather than identical, noise cancellation, for the reason 
of further reducing the costs. In this case, optimal solutions that still guarantee the exact cancellation are 
likely to provide good initial guesses for subsequent optimization in the approximate framework. 

As indicated previously, in the current paper we focus primarily on the quadratic optimization criteria. 
We have looked into the most natural criterion of this type, namely, the L 2 norm of the control sources g(x) 
of (2.1) or gffi of (2.11), see Section 3.2. Clearly, this cost function depends only on the controls themselves. 
The advantage of minimizing the controls in the sense of L 2 : 



is that the minimum can be easily computed , see Section 3.2. The search space for minimization (3.1) includes 
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all the appropriate auxiliary functions w(x), on which g(x) depends. The disadvantage of using this criterion 
is that the quantity ||p||2 does not have a clear physical interpretation. Nonetheless, motivated primarily 
by the ease of the numerical approach to minimization, we do provide in Section 3.2 a comprehensive set of 
computed optimal solutions for active controls in the sense of the least squares (i.e., L2). We also compare 
these discrete results with the “semi-analytic” L2 -optimal solutions obtained for simple circular shapes using 
the spectral methodology developed in our previous paper [4]. 

Note, an alternative to minimization in the sense of L 2 (3.1) may be minimization in the sense of L\: 

Mil = J > min . 

supp 9 

We have thoroughly studied this problem in our recent paper [5]. In particular, we have shown that the 
L\ minimization is equivalent to minimizing the overall absolute acoustic source strength , see [6,7], of the 
control sources g{x). This clear physical interpretation constitutes an advantage of using the L\ norm of the 
control sources as a cost function for optimization (besides that it also depends only on the control sources 
g( x) themselves). On the other hand, the corresponding optimization problem has proven difficult to solve 
numerically, see [5]. We briefly describe the L\ results of [5] in Section 3.1 for the purpose of comparison. 

In the discrete framework, the minimization problem for the control sources can be formulated as 
follows: 

min, 

and the L\ minimization problem can be formulated as follows: 

lk (h) ||i= H VMSm’l — unm, 

mfrM - 

where V m accounts for the cell area and again, the search space includes all the appropriate auxiliary grid 
functions w^ h \ through which g ^ is defined, see formula (2.11). 

Either of the two foregoing discrete minimization problems can also be rewritten using matrices. The 
finite-difference operator L W can obviously be interpreted as a matrix with N columns and M rows, where 
N is the number of nodes n = (s, j) of the grid N _ such that pj < ii, i.e., j < J, and M is the number of 
nodes m = (s,j) of the grid Mr such that pj < it, i.e., j < J - 1. Let to be the vector of N components 
w = w[ k j , n £ N“ and j < J, arranged so that 

w = [tf> 7 , too, to ,J-1, to.,j] T , (3.2) 

where to 7 contains for which n € 7, w.,j and to.,j_i correspond to the outermost and second to last 
circles of the polar grid, respectively, as in formula (2.15), and too contains all the remaining components of 
to “in-between” 7 and the outer boundary. In accordance with (3.2), the matrix X^) can be decomposed 
into four sub-matrices: 

L {h) = [A, B, C, D] (3.3) 

that all have the same number of rows M, A has as many columns as there are nodes in 7 (we denote this 
number I7I), C and D each has L columns, and B has A7 — (7! — 2 L columns. 



13 



Using formulae (3.2) and (3.3) and introducing a generic notation || * || for either || ■ ||i or || • || 2 , we have: 

|| V(Aw 7 + Bwo + Cw. j- 1 + Dw. y j)\\ — > min, (3.4) 

where V is an M x M diagonal matrix with the entries given by the corresponding cell areas V 7 m . The 
vector w in the optimization formulation (3.4) is, in fact, subject to a number of equality-type constraints 
that come from the interface conditions (2.12) and ABCs (2.15). More precisely, the first sub- vector iu 7 in 
(3.2) is known and fixed because of (2.12) and we can rewrite (2.12) as w 1 = where u 7 is given. The last 

sub- vector w.j in (3.2) is a function of according to (2.15). Therefore, we can conclude that only tub 

and w. j-i contain free variables that provide the search space for optimization, and as such rewrite (3.4) as 

min || V(Bwo + (C + DT)w. }J -i + ^w 7 )|| = min||fk -/||, (3.5) 

U70,tiJ. ( j— j Z 

where E = V[B, C + DT] is an M x (TV — |^y | - L) given matrix, z = [«*), w. t j-i] T is an (N — |^y| — L)- 
dimensional vector of unknowns, and / = - VAw y is an M-dimensional known vector of the right-hand side. 
Minimization problem (3.5) is, in fact, a problem of finding a weak solution of an overdetermined complex 
linear system Ez — /. Note, the quantities involved are complex because we are dealing with the traveling 
waves in the framework of the Helmholtz equation (see Section 1). 

3.1. Results of Optimization in the Sense of L\. To provide a “reference point” for comparison, 
we outline here the findings of our recent work [5]. Because of the complex- valued quantities involved, the 
computation of the weak solution of problem (3.5) in the sense of L\ reduces to solving a non-linear and non- 
smooth problem of constrained optimization over a large set of cones. This problem presents a substantial 
challenge even for the most sophisticated state-of-the-art approaches to numerical optimization, namely, 
those based on interior point methods [8, 16]. The difficulties are further exacerbated by the large dimension 
of the grid, because on one hand, the number of conical constraints that one obtains when solving (3.5) in 
the sense of L\ is the same as the number of nodes in the grid M, which can be quite large even in two space 
dimensions, and on the other hand, the typical maximum number of constraints that the aforementioned 
state-of-the art methods can handle is only on the order of hundreds. 

In spite of the difficulties, we have been able to compute several two-dimensional solutions for simple 
test cases. Our best results were obtained with the software package SeDuMi by J. F. Sturm. 4 This is a 
numerical algorithm for optimization over cones [13], it employs the ideas of interior-point methods, and the 
self-dual embedding technique of [19], see also [9]. The algorithm allows for complex- valued entries, which is 
very important in our framework, and also for quasi-convex quadratic and positive semi-definite constraints. 

All numerical experiments that we have conducted, see [5], indicate a very consistent behavior of the L\ 
optimal solution for the control sources. It happens to be the discrete layer of surface monopoles tf^onopoie 
described in the concluding part of Section 2. Recall, this solution is obtained by applying formula (2.11) to 
the auxiliary function w W defined by (2.16a), (2.16b). Figure 3.1 reproduces the results of one particular 
computation from [5]; the protected region ft for this case was a unit disk, and the grid dimension was chosen 
32 x 7 (see Section 3.2.1 for the actual grid definition). 

Motivated by the consistent observations of the L\ optimal solution being the same as surface monopoles 
0monopoie> we h ave a ^ so b een a ^ e to prove in [5] that surface monopoles indeed provide a global minimum for 

4 http : //f ewcal . kub . nl/ Sturm/ oof tvare/sedumi . html 
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(a) L\ optima] solution 


(b) Surface monopoles S^onopL 


Fig. 3.1. Magnitude of the control sources for ft = {* € R 2 ||x| < 1}, A: = 0.5, / = S(x — *i), Xi — (5,0). 


the control sources in the sense of Li for both the discrete and continuous formulation of the problem. The 
proof of [5], however, covers only the one-dimensional case. Even though we have not yet been able to prove 
a similar result for the general multi-dimensional case, we still believe that it is true, because a combination 
of the two-dimensional numerical evidence and a one-dimensional accurate proof cannot, in our opinion, be 
a mere coincidence. Therefore, we formulated this result in [5] in the form of a conjecture that we reproduce 
below. Let us remind, that according to (2.3) the continuous surface monopole controls are given by 

- - [Si - - ^]/ <r| - -wLr ■ <<r) ' (3 - 6) 

where iu(x) is a solution to the exterior Dirichlet problem (2.4). Then, we have 

CONJECTURE 3.1. Let a complex-valued function w = w(x) be defined on Hi = and let it be 

sufficiently smooth so that the operator L of (1.1) can be applied to w(x) on its entire domain in the classical 
sense , and the result Lw be locally absolutely integrable. Let , in addition u>(x) satisfy the interface conditions 
(2.2), where u = u(x) is a given field to be controlled , and the appropriate Sommerfeld radiation boundary 
conditions at infinity , (1.2a) or (1.2b). Then , the greatest lower bound for the Li norms of all the control 
sources g(x) obtained with such auxiliary functions w(x) using formula (2.1), is given by the L\ norm on T 
of the magnitude of surface monopoles (3.6): 

inf [ \g(x)\dx = f \u{x)\ds. (3.7) 

«>(*) JOa Jr 

Formula (3.7) can obviously be rewritten as 

inf ||ff(*)||i,o 1 = IMkr* ( 3 - 8 ) 

w{x) 

As has been mentioned, the discrete prototype of formula (3.8) that reads as follows 
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was conjectured in [5] on the basis of the two-dimensional experimental observations. We emphasize that in 
the discrete case the L i -optimal solution ^onopoie belongs to the same class of functions as all the discrete 
volumetric controls of (2.11), whereas in the continuous case the optimum on the class of volumetric 
controls g(x ) of (2.1) actually takes us out of this class to the singular layer (3.6) on the interface T. 

3.2. Discrete Optimization in the Sense of L 2 . The L 2 minimization problem for the volumetric 
control sources is solved hereafter completely on the discrete level. In other words, for every particular setup 
we are finding the minimum (3.5) or, equivalently, computing a weak solution of the overdetermined system 
of linear equations Ez = /, in the sense of the least squares. The resulting optima do not reduce to any 
clearly identifiable special cases, like the layer of surface monopoles that appeared in the previously analyzed 
context of L\. They are not assigned any particular physical meaning either, we present them below in order 
to demonstrate that the L 2 optima are distinctly different from the L x optima obtained in [5], and that they 
can be easily computed numerically, including some cases that involve rather sophisticated geometry. In the 
simple case when the protected region fi is a disk, we also conduct a grid convergence study in order to 
validate the results of the discrete L 2 minimization against the analytic reference solutions computed with 
high accuracy using the spectral methodology that was first proposed in [4]. 

Proposition 3.1. The matrix E = V[B, C + DT], see formulae (3.3), (3.5), has full column rank. 

Proof. The justification of Proposition 3.1 will be based on a natural solvability assumption for the 
system of finite-difference equations that we are using. First, let us introduce more detailed partitions of w 
and L^' 1 instead of (3.2) and (3.3), respectively: 

w = [tu 7 +, u> 7 - , Wq, tO.'J-i, w.'j] T , 

(3 9 ) 

L (h) = {A + ,A~,B,C,D]. 

The matrices A+ and A of (3.9) together give A of (3.3); tw 7 + and A+ correspond to the innermost 
“half” of the grid boundary 7 +, and u> 7 - and A~ correspond to the outermost “half” of the grid boundary 
7 “ (see formula (2.10) and Figure 2 . 1 ). Next, consider an auxiliary exterior Dirichlet problem for the 
finite-difference equation = 0 [see formula (2.9)] with the boundary data specified at 7 +. As 

before, the problem is supposed to be truncated at the external artificial boundary p — pj by means of 
the ABC (2.15). This problem is a discrete counterpart of the continuous exterior Dirichlet problem for 
the Helmholtz equation with the boundary data given at T and ABCs (2.7) specified at p = R, The 
continuous problem is uniquely solvable because it is equivalent to the genuine infinite-domain exterior 
Dirichlet problem with the Sommerfeld boundary conditions (1.2b) set at infinity. Even though we do not 
prove it, it is certainly reasonable to assume that the corresponding discrete problem based on a standard 
central- difference scheme (2.9) and ABC (2.15) is uniquely solvable as well . 5 The latter assumption implies 
that the square M x M matrix C + DT], see formula (3.9), is non-singular. Consequently, the 

matrix G — V[A , B , C DT] is also non-singular, because V is an M x M diagonal matrix with non- 
zero diagonal entries V m . Finally, we notice that the matrix E = V[B , C + DT ], see formulae (3.2) and 
(3.3), is obtained by removing the first | 7 ~| columns of the previous matrix G. Therefore, the columns of 
E are linearly independent. □ 

5 A proof of this fact would involve showing that relations (2.14) axe “sufficiently close” to guaranteeing the precise mode 
selection in the discrete case so that to avoid the resonances. This task is beyond the scope of the current paper. A comprehensive 
study of solvability and well-posedness of one-dimensional discrete boundary value problems can be found in [11]. 
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An obvious key implication of Proposition 3.1 is that the minimization problem (3.4), or equivalently 
(3.5), can be solved in the sense of L 2 (least squares) using a standard QR-based approach, i.e., without 
employing the Moore-Penrose type arguments. We use the MATLAB function LSQLIN for solving the least 
squares minimization problems hereafter. This function also allows one to do constrained minimization, the 
capability that we employ in Section 3.2.2. 

3.2.1. Comparison with the Analytic Solution. In our previous work [4] we have developed a 
methodology of spectral type that allowed us to construct the continuous Z/ 2 -optimal volumetric controls 
for a particular geometry, namely, controls supported on annular domains. We employed the separation of 
variables and expressed the exact optimum as an infinite Fourier series in the circumferential direction whose 
coefficients were certain combinations of Bessel functions, see [4, formula(5.21)]. This series obviously had 
to be truncated at a certain maximum number of harmonics for the purpose of numerical evaluation. On 
smooth solutions, this method obviously provides for a spectral convergence. For the current purpose of 
validating the finite-difference algorithm we will use the spectral solution of [4] as a reference solution in the 
grid convergence tests. 

Let the protected region be a disk of radius r centered at the origin: ft = {(p, 0)\p < r}, and let the 
controls be supported on the annulus Hi = {(p,0)|r < p < R}. We introduce a simple conformal polar grid, 
which is uniform in the circumferential direction and stretched in the radial direction so that the cell aspect 
ratio is equal to one: 

M={(p j ,9,)\p j =e^ e ,j = 0,... ,J- 1; 0 g = sA0,s = 0, . . . ,L - 1; AO = 2i r/L}, 

N = {{pj,O s )\ pj = e jA9 , j = —1,0,... , J; 0, = sA0,s = 0,... ,L- 1; AO = 2tt/L}. 

We, of course, assume that the area covered by the grid N of (3.10) is larger than fli, i.e., p_i < r < R < pj. 
The Helmholtz operator can be easily approximated on the grid (3.10) with the second order of accuracy 
using the same five-node stencil as shown in Figure 2.1. This approximation involves only minor changes 
compared to the approximation (2.9) that works on uniform grids, and we refer the reader to our paper [10] 
for detail. The discrete ABCs in the form (2.14) or (2.15) do not change, except that A p needs to be replaced 
by A pj = pj — pj-i . Let us also note that we do not consider the grid (3.10) inside the domain fl because 
we introduce it only for the purpose of obtaining the control sources on fli. If, however, we were to actually 
compute the output of the controls inside the protected region, we would have had to extend the grid (3.10) 
all the way into fl, which can obviously be done using a variety of strategies. As indicated by the previous 
analysis [4,12,17,18], as long as the discrete controls are constructed according to formulae (2.11), (2.12), 
their output on the grid inside fl will identically cancel out the unwanted acoustic component see 

Section 2.2. In all the cases that we analyze hereafter, we have p_ i < r < po, so that according to (2.10) 
the grid subsets are defined as 

M+ ={(Pi,0.)|i = -l}, M- ={(pM\ 0< j< J-l}, 

ft = N- ={(P,A)| - l<j< J), (3.11) 

7 = {(Pji^s) \j = —1)0}, 7 + = {{PjJ»)\i = -1}> 7" = {(PjA)|j =0}, 

where always a = 0, ... ,L — 1. In so doing, the dimension of the matrix L^ h \ see (3.3), is M x .V = 
(L • J) x (L • ( J + 2)), the dimension of A, which corresponds to the variables on 7, is M x 2 • L ~ (L « J) x 2 • L, 
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the dimension of B is M x (N - 4 L) = (L • J) x (L • ( J — 2)), and the dimension of either C or D is 
M x L = (L • J) x L. 

We test the convergence of the discrete scheme for the wavenumber k — 1.0 and the excitation (i.e., 
the acoustic field u ^ that drives the control system) taken in the analytic form of a shifted fundamental 
solution of the Helmholtz operator, as if it were generated by the point source f~ = 5{x — X \ ), where 
X\ = (p cos 0, p sin 0) = (5,0). We reemphasize that our approach does not require the explicit knowledge of 
the exterior sources of noise. We only need this function u ^ as a sample field to be used as given data in 
formula (2.12). 

We employ a sequence of seven grids: Lx J = 32 x 3, 48 x 4, 64 x 5, 96 x 7, 128 x 9, 192 x 13, and 256 x 17, 
so that for all the grids the value of pj-\ is the same: = const « 1.481; according to (3.10) we also have 

po — 1. For the first series of convergence tests we assume that the boundaries p = r and p — R of the region 
fli, on which the continuous controls are to be supported, are located exactly at the conformal midpoint of 
the first and last cell of the radial grid N of (3.10), respectively, i.e., r = e~ l / 2Ae and R = The 

results of these tests are summarized in Table 3.1, which shows the L 2 norm of the relative error between 
the optimal continuous and discrete controls: arg[min ||p S pect(^) II 2 ] and arg[min \\g^ lli^l . 

L w{x) L u>(M 1 z 


Table 3.1 

Grid convergence for: k = 1 , /” = S(x — Xi), X\ — ( 5 , 0 ), r — R = 


Grid 

32 x 3 

48 x 4 

64 x 5 

96 x 7 

128 x 9 

192 x 13 

256 x 17 

||Error|| 2 

0.013722 

0.0061417 

0.0034693 

0.0015491 

0.00087349 

0.00038921 

0.00021922 


The data in Table 3.1 clearly indicate the second order of grid convergence for the discrete optimal 
controls g^ h \ It is important to emphasize, though, that the geometry of fii was chosen grid dependent 
(boundaries p — r and p — R were located at cell midpoints), which essentially means that for each subsequent 
grid in Table 3.1 the optimum was computed on a somewhat different (smaller) domain. It is quite obvious 
that in general the optimal solution will depend on the region on which the optimization is performed, and 
we cannot expect the optimum computed on a subdomain to coincide with the corresponding fragment of 
the optimum computed on the entire domain. However, the decrease of the error with the refinement of the 
grid observed in Table 3.1 shall still be interpreted as convergence. Indeed, had we continued refining the 
grid further, all the domains Hi = {r < p < i?} themselves would converge to one and the same annular 
region with the inner radius r = po = l (Q is & unit disk) and outer radius R = pj_i, which was chosen grid 
independent. 

On the other hand, the quadratic rate of convergence suggested by Table 3.1 appears a rather fragile 
phenomenon determined by the particular choice of the geometry. For other choices, the convergence may 
be slower. In Table 3.2, we present the results that correspond to the same inner boundary r = e~ 1 / 2A6> , 
and the outer boundary located at either one quarter point or three quarters point of the outermost cell: 
R = e ( J ~ 3/4) Afl or R = One can easily see that the convergence in Table 3.2 is only linear. 

At the moment, we do not have a detailed explanation of the grid convergence properties for g W that 
we have observed, see Tables 3.1 and 3.2. It is important to realize, however, that what we evaluate is, in 
fact, convergence of the residual rather than that of the solution. Indeed, the solution of the optimization 
problem (3.4) or (3.5) per se is a particular grid function w ^ that delivers minimum to the selected function 
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Table 3.2 

Relative L 2 error for: k = 1, = S(x — xi), x\ = (5, 0), r = , R -- e^* 7-3 / 4 )^ and R = e^ J 1 /4)A0 , 


Grid 

32 x 3 

48 x 4 

64 x 5 

96 x 7 

128 x 9 

192 x 13 

256 x 17 

J — 3/4 

0.10283 

0.074135 

0.058274 

0.040958 

0.031618 

0.021733 

0.016562 


0.096894 

0.071804 

0.057069 

0.040476 

0.031362 

0.021627 

0.016504 


of merit, namely, the L 2 norm of the residual of the discrete Helmholtz operator applied to this What 
motivates our primary interest toward the residual is obviously the fact that it has a physical meaning of the 
distributed active control sources g( h \ see formula (2.11). However, from the numerical analysis standpoint 
it is known that grid convergence of the residuals is, generally speaking, not guaranteed even if the solution 
itself does converge. Moreover, even though the optimization formulation that we have introduced in the 
beginning of Section 3 is fairly conventional, in the PDEs’ perspective neither the continuous generating 
function w(x) nor its discrete counterpart w (*> at the optimum can be interpreted as a solution to any 
traditional boundary- value problem, for which the existence and regularity results are available. As such, no 
standard theoretical approaches to analyzing grid convergence will directly apply here, and we shall rather 
regard the foregoing results as experimental findings. 

Let us point out that in the context of noise cancellation on the domain ft, the issue of grid convergence 
of the discrete control sources gW may, in some sense, be considered as the one of secondary importance. 
Indeed, the output of the controls g W always eliminates the unwanted noise on ft [more precisely, on the 
grid N + , see formula (2.10)] no matter what particular solution from the general class (2.11), (2.12) is used. 
Moreover, this output on N + can be interpreted as a discrete generalized potential of Calderon’s type, which 
will always converge to its continuous counterpart with the rate prescribed by the approximation order of 
the scheme, again, irrespective of what particular , n E N”, and \ m E M” , are taken to generate 

the potential on every given grid, see the discussion on page 11 of this paper and references [4, 12] for 
detail. As such, one need not be overly concerned with the rate of convergence for the discrete optimal 
control sources as any of those will do the cancellation job equally well in any event. The question of grid 
convergence for g^\ however, may be of a considerable independent interest, from both the theoretical and 
experimental standpoint. It may certainly be worth looking into in the future, even though we expect that 
neither theoretical nor systematic experimented analysis will be straightforward, especially in the case of 
general geometries. The focus of the current paper, however, is not so much on the study of grid convergence 
for g^ h \ but rather on building and testing the discrete quadratic optimization algorithm and applying it to a 
number of cases, including those for which little is known regarding the analytic solution (see Section 3.2.2). 

What we also want to emphasize in the current paper is that the L 2 optimal solutions for active controls 
differ very substantially from the L\ optimal solutions obtained previously in [5]. This is, of course, natural 
to expect, but it is also interesting to visualize and actually observe the corresponding differences. As such, 
we proceed with conducting the least squares minimization for the same setup that was earlier analyzed in 
the sense of L\ in our paper [5] (the L\ results are reproduced in Figure 3.1). The grid dimensions were 
L = 32 and J — 7; and the wavenumber k in the Helmholtz equation (1.1) was chosen k = 0.5. The 
excitation was again produced by the point source f~ — 6(x - X\) y where X\ = (5,0). In Figure 3.2(a) 
we show the magnitude of the L 2 -optimal active controls on the 32 x 7 grid. This solution indeed differs 
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drastically from the Li-optimal controls that are shown in Figure 3.1. Unlike the L\ optimum, i.e., the layer 
of surface monopoles, the L 2 optimal solution tends to be distributed over the entire annular region on which 
the control sources are supported, obviously favoring the direction toward the noise source. We have also 
obtained the L >2 optimal controls for the same case but on a twice as fine grid of dimension 64 x 13; they are 
shown in Figure 3.2(b). The plots in Figures 3.2(a) and 3.2(b) look very much alike, as expected. 


L2 optimal control offort jg| 


Analytic L2 optimal oonboi aflort |q{ 



(a) L = 32, J = 7 (b) L = 64, J - 13 

Fig. 3.2. Magnitude of the I- 2 -optimal control sources for fl = {i E R 2 1 |x| < 1}, k = 0.5, f~ = 6(x — xi), xi = (5,0). 


It is also interesting to observe how the qualitative behavior of the optimal solution changes when the 
parameters that define the problem change. A key parameter is the wavenumber k. Previously, we have 
analyzed the cases of relatively long waves compared to the size (i.e., diameter) of the protected region fl. 
Let us now take k — 7r, then there will be exactly one full wavelength across the diameter. We compute this 
case on the grid 128 x 9 so that 1 = po < p < pj-i « 1.481. In Figure 3.3, we present the distribution of 
optimal controls g ^ for the case of the long waves, k = 0.5 [Figure 3.3(a)], and for the case of the wavelength 
comparable to the domain size, k = n [Figure 3.3(b)]. One can clearly see that the solution that corresponds 
to shorter waves is more oscillatory. 

3.2.2. Constrained Optimization in the Sense of L2* The purpose of formulating and solving the 
L 2 optimization problems that involve constraints was to simulate not simply a more sophisticated geometry 
but also a more realistic one. For example, if we interpret the previously considered protected region — a 
unit disk — as a section of the aircraft fuselage, then we can also introduce portholes, i.e., windows, that 
shall be interpreted as designated areas, in which no control sources can be applied. Optimization problem 
(3.5) in this case needs to be modified. Instead of simply finding a weak solution of Ez = / in the sense 
of the least squares, we will now have to impose additional constraints, i.e., require that for those nodes of 
the grid that happen to be inside the aforementioned designated areas, the corresponding equations be 
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L2 op*m*l control •ftort |g| 


12 optimal control *flort lot 




(a) k = 0.5 

Fig. 3.3. Magnitude of the L 2 -optimal control sources for S7 = {x € K 2 ||®| < 1}, f~ = 5{x - ®i), x\ = (5,0), 128 x 9 grid. 


enforced exactly. This leads to the problem 

min H-Ez - /H 2 subject to E c z = J, z , (3-12) 

2 

where E c is the sub-matrix of E (i.e., the appropriate set of rows), and f c is the respective sub- vector of /, 
that correspond to the constrained nodes. 

For simulations, we have introduced two symmetrically located portholes in the fuselage: 5° < 8 < 30° 
and 150° < 6 < 175°. The resulting problem (3.12) was solved by a standard methodology for the least 
squares minimization with equality constraints. It requires that the constraints be linearly independent and 
basically results in reducing the dimension of the remaining search space accordingly. In our computations, 
we have used the procedure LSQLIN available in MATLAB. 

The case that we have actually analyzed in the context of the constrained L 2 optimization, was again 
one of those that we have studied previously in the L\ framework, see [5], but obviously with no constraints. 
For this case, the excitation is provided by a pair of external sources: /” = S(x - X\) + - afc), where 

x\ = (5,0) and x^ — (1,2), the wavenumber k = 0.9, and the original grid has the dimension 48 x 9. In 
Figure 3.4(a), we show the constrained L 2 optimal solution for this original grid, and in Figure 3.4(b) we 
show the solution for the twice as fine grid 96 x 17. We emphasize the presence of the large spikes in the 
control effort next to the boundaries of the window on the right, which is natural to expect. We should also 
point out at some apparent discrepancies between the control field on Figure 3.4(a) and that on Figure 3.4(b) 
in the region near this window. Qualitatively, these discrepancies are easily explained once we realize that a 
given window, which is defined as a particular range of 0, does not have to be exactly the same on different 
grids because of the finite size A0, and a finer grid simply provides for a “sharper” definition of the window in 
the discrete sense. On the other hand, quantitatively we, of course, cannot claim that the same convergence 
results as we have obtained previously in the case with no constraints, see Section 3.2.1, will hold in the 
presence of the constraints as well. Moreover, in the constrained case one should generally expect less 
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(a) L = 48, J - 9 


(b) L = 96, J = 17 


Fig. 3.4. Magnitude o/ £/ie L? -optimal control sources for ft = {x 6 R 2 ||*| < 1}, k = 0.9, excitation: / = <5(ac — xi) + 
<f(x — Z2)> asi = (5,0), ®2 = (1,2), and window constraints: 5° < 6 < 30° and 150° <6 < 175°. 


regularity from the corresponding continuous solution than in the previously addressed unconstrained cases. 
Therefore, the results of the L 2 constrained minimization outlined in this Section should only be regarded 
as implementation examples of a previously tested numerical algorithm for more elaborate settings. 


Table 3.3 

Comparison of the computed L 2 -optimal solutions with surface monopoles. 


Grid 

min g {h) (h) 

tu<*> 2 

Constrained min 

(h,surf) I (*) 
^monopolef 2 

48 x 9 

0.41855 

0.54013 

1.2983 

96 x 17 

0.43485 

0.56175 

1.8315 


It is also interesting to compare the actual norms of the solutions that we have obtained. They are 
presented in Table 3.3, which also contains the L 2 norms of surface monopoles that are optimal in the sense 
of Li, see [5]. From Table 3.3 we see that the L 2 norm at the minimum is considerably larger for the 
constrained case compared to the unconstrained case. As concerns the L 2 norm of the L\ -optimum, it is 
three times larger in this case than the unconstrained L 2 minimum. We should also mention that the finer the 
grid, the larger the L 2 norm of Pmonopoie see Table 3.3. This is, in fact, a natural consequence of the scaling 
that we have adopted in [5]. Indeed, as indicated in [5], the actual magnitude of ^onopoie i ncreases when 
the grid is refined, because the corresponding continuous limit is a single layer on the interface. The latter is 
a singular distribution, which is obviously not integrable by itself, and even less so with square. At the same 
time, it turns out that the discrete two-dimensional L\ norm of surface monopoles ||Pmonopoieili m- does not 
change with the change of the grid size. This essentially implies that the magnitude of P^onopoie sca ^ es as 
0(h~ l ) and as such, the L 2 norm ll^onopoiel^ m- * s su PP ose d to scale as 0(h ~ 1 / 2 ). This is corroborated 
by the data in the last column of Table 3.3. 
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4. Discussion. We have developed and implemented a computational algorithm for optimizing the 
sources of active control of sound in the sense of L 2 . For simple cases, we have been able to validate our 
numerical results against the analytic solution. We have also seen that the L 2 optimal controls are distinctly 
different from the Li optimal controls obtained previously. For the case of a somewhat more realistic 
geometry, the corresponding optimization formulation involves constraints of equality type. Our algorithm 
allows us to analyze the constrained L 2 optimization problems as well. 

In general, we should mention that there is a multitude of different acceptable optimization criteria for 
active control of sound. For example, the advantage of L 2 is its clear physical interpretation as minimization 
of the overall absolute acoustic source strength. L 2 does not have such a transparent physical meaning, 
but is easier to compute numerically. In the forthcoming paper 6 , we will report the results of the power 
optimization. It turns out that the corresponding analysis necessarily involves interaction between the 
sources of sound and the surrounding acoustic field, which is not the case for either L\ or L 2 . Even though 
it may first seem counterintuitive, one can build a control system that would require no power input for 
operation and would even produce a net power gain while providing the exact noise cancellation. Of course, 
other functions of merit, besides the aforementioned three, can be employed as well. Some may come from 
the engineering limitations, others wall be just a matter of personal preference. Questions related to the 
optimization of active controls of sound using different criteria will be studied in the future. 
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